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Large-eddy simulations with wall models 

By W. Cabot 

1. Motivation and objectives 

The near- wall viscous and buffer regions of wall-bounded flows generally require 
a large expenditure of computational resources to be resolved adequately, even in 
large-eddy simulation (LES). Often as much as 50% of the grid points in a com- 
putational domain are devoted to these regions. The dense grids that this implies 
also generally require small time steps for numerical stability and/or accuracy. It 
is commonly assumed that the inner wall layers are near equilibrium, so that the 
standard logarithmic law can be applied as the boundary condition for the wall 
stress well away from the wall, for example, in the logarithmic region, obviating 
the need to expend large amounts of grid points and computational time in this 
region. This approach is commonly employed in LES of planetary boundary layers 
(e.g., Mason, 1989; Schmidt &: Schumann, 1989), and it has also been used for some 
simple engineering flows (e.g., Piomelli et a/., 1989; Arnal & Friedrich, 1993). 

In order to calculate accurately a wall-bounded flow with coarse wall resolution, 
one requires the wall stress as a boundary condition. The incompressible Navier- 
Stokes equation is 


— = -Vp + Vr, t = — uu + i/Vu , (1) 

in which u is the velocity, p is the pressure, r is the stress, and v is the molecular 
viscosity. In a simulation with an unresolved wall, the wall-normal (y) derivative of 
the stress for tangential ( x,z ) velocity components, 

Ty + "f^f) • i = 1 ’ 3 ' < 2 > 

cannot be accurately calculated by applying the usual no-slip condition, u = 0, 
instead requiring the specification of the wall stress 

, * = 1,3. (3) 

y = o 

Thus, an adequate model of based on outer flow quantities is desired. Asymp- 
totic matching of inner and outer regions in steady, ensemble-averaged, equilibrium 
flow yields the log-law relation between wall stress and outer mean velocity. How- 
ever, for the purposes of LES, wall stress models are needed with some degree of 
time and space dependence. Because the near-wall layer is typically very thin with 
respect to horizontal scales, boundary layer assumptions may be valid, perhaps even 


Ti2w = V 


dui 

dy 
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on horizontal grid scales used in LES, and it may be possible to use simpler bound- 
ary layer equations to model the near- wall region and at the same time retain more 
flexibility in handling flows with widely varying pressure gradients. 

The goal of this work is to determine the extent to which equilibrium and bound- 
ary layer assumptions are valid in the near-wall regions, to develop models for the 
inner layer based on such assumptions, and to test these modeling ideas in some 
relatively simple flows with different pressure gradients, such as channel flow and 
flow over a backward-facing step. Ultimately, models that perform adequately in 
these situations will be applied to more complex flow configurations, such as an 
airfoil. 


2. Accomplishments 

An examination of momentum balance at different horizontal scales, and corre- 
lations between the measured wall stress and some outer flow quantities, have been 
performed from a direct numerical simulation (DNS) database for channel flow. Be- 
cause wall stresses need to be predicted in flows with different pressure gradients 
and in separated flow, models based on the log law and boundary layer equations 
have been tested both in channel and backward-facing step flows. 

2.1 Momentum balance in channel flow 

Near-wall data has been examined from a channel flow DNS (J. Kim, personal 
communication; Kim, Moin Sz Moser, 1987) with a friction Reynolds number Re T = 
395 ( Re r = u t 8/v , where 6 is the channel half-width, u T = \vdll /dy\ l l 2 is the 
friction speed, and U is the mean streamwise velocity). Horizontal averages of 
flow quantities were taken over different scales, from the scale of the entire plane 
down to scales comparable to expected LES resolutions (a factor of 16 smaller in 
each direction, or x A « 160 x 80 in wall units scaled by v/u r ). The 

streamwise momentum balance was constructed by integration over volumes with 
these horizontal dimensions from the wall to a height y + ^ 80: 


du du 2 duw dp ' \ 

dt + !h + ~dT + / 


-f x +(*•<■*■)- 


duv\ 

W/ ’ 


( 4 ) 


where (• * •) denotes a volume average, and dP/dx is the mean pressure gradient. 
The results shows that the advection and fluctuating pressure gradient terms on 
the left-hand side of (4), while small compared to the other terms when averaged 
over the entire plane, are more than an order of magnitude larger at LES scales. 
This suggests that momentum balance is dominated by a nearly inviscid balance 
between advection and pressure gradients at LES scales, casting doubt on the local 
validity of models, such as the log law, based on a balance between terms on the 
right-hand side of (4) (J. Jimenez, personal communication). 

Correlations between the wall stress T\2 W and the mean streamwise velocity at 
y + ss 40 are small but significant (50% at LES scales). Figure 1 shows a scatter 
plot of the deviation from the mean of actual wall stress versus that predicted from 
a logarithmic law with the (nearly zero) mean pressure gradient in the channel at 
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FIGURE 1. Deviation from the mean of actual wall stress from DNS channel flow 
data (Re r = 395) compared with that predicted from the log law model applied at 
y + « 40. The flow is averaged horizontally on typical LES scales (A:r+ x A « 
160 x 80). The linear diagonal denotes a perfect local correlation. 

ps 40. There is a noticeable linear correlation for values of wall stress near the 
mean, with larger deviations in high-stress regions. (The nature of the high-stress 
events has yet to be explored.) On the other hand, the correlation of wall stress to 
the large, instantaneous, fluctuating pressure gradients is found to be practically nil 
(only a few percent). Corresponding analyses need to be performed with DNS and 
LES databases for flow over a backward-facing step (Le h Moin, 1993; Akselvoll & 
Moin, 1995), which contain a large adverse pressure gradient and separated flow. 

2.2 Boundary layer wall models in channel flow 

Wall models have been tested in a second-order, central finite difference (FD2) 
channel code on a staggered mesh with a third-order Runge-Kutta (RK3) time 
advancement (Akselvoll & Moin, 1995), in which the wall stress boundary conditions 
are easily implemented. These wall models have been based on the Johnson-King 
(1985) boundary layer model, which is fairly simple and has had good success in 
Reynolds- averaged Navier-Stokes (RANS) models of separated flow (Menter, 1991). 

A channel flow with a target Re r = 1030 was simulated using an outer mesh with 
the near-wall points for horizontal velocity placed at a matching height = 32 
or 64. Embedded in the outer mesh is a fine sublayer mesh from the wall to the 
matching height. The outer mesh technically extends to the walls, but only the 
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v = 0 and dp/dy = 0 boundary conditions are used. Both outer and inner meshes 
are usually stretched with a hyperbolic tangent mapping. The outer mesh uses 33 
wall-normal nodes, and the sublayers uses 21 nodes at each wall. The horizontal 
domain size is Ax x A z = 2ir6 x irS. Initially, a horizontal mesh of 32 x 32 was 
used for both outer and sublayer regions, but it was found that much better results 
were obtained with a finer 64 x 64 mesh for the outer region; on the other hand, 
the mean velocity and rms statistics were found to be insensitive to whether the 
sublayer mesh was 32 x 32 or 64 x 64 (Ax + x A z 4- « 200 x 100 or 100 x 50). It 
was also found that results from the FD2-RK3 code were sensitive to the time step 
for convective CFL numbers exceeding about 0.5, perhaps due to inaccuracies in 
implicit terms (cf. Choi & Moin, 1994). In the results presented here, the convective 
CFL number was kept around 0.6. 

Model JK0. The lowest level model for the wall stress is obtained at each hori- 
zontal position in the near-wall sublayer (independent of other horizontal locations) 
from the solution of the ordinary differential equation 


d dU, dP 

-y{v + Vt)— = — , i = 1,3 , 
ay ay axi 


(5) 


where Z7, are the horizontal velocity components in the sublayer, dP/dxi is the 
constant mean pressure gradient, and 

v t = Ku s y w D 2 , D = 1 — exp(-Udyw/Au) , (6) 


resembles the eddy viscosity in the Johnson-King (JK) model for the inner regions. 
Here, though, the scale speeds u 3 and Ud are replaced by the friction speed u T \ y w 
is the distance from the wall, k is the von Karman constant, and A is a damping- 
function constant taken to be 19, which gives the best fit to the standard log law in 
this case (lower values were used by Johnson & King and Menter). The boundary 
conditions for (5) are Ui — 0 at the walls and 17, equal to the horizontal velocity in 
the outer mesh at the first grid point above the wall. The wall-normal derivative 
of Ui at the wall yields the wall stress r^w used in the outer flow. Eq. (5) is solved 
by using the same FD2 discretization used in the main code and performing an 
inversion of the resulting tridiagonal matrix. The solution of (5) is just a smooth 
blend of the viscous and logarithmic functions and, for the channel, is generally 
equivalent an instantaneous log law. Because one can consider expressions like (5) 
to be valid only in some average sense, both in space and time, a running time- 
average of the matching velocity over about an eddy turnover time is employed. 

Model JKOa . The next level of model tests the influence of large advective and 
instantaneous pressure gradient terms: 


d, , dUi dU t „ /rrTT , dp 

-«r 


* = 1,3 , 


(7) 


where Ui are the horizontal velocity components in the sublayer, as in (5). The 
eddy viscosity is given by (6) and, as in the JK0 model above, uses u a = Ud — u r . 
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The solutions of (7) at different horizontal locations are now coupled through the 
divergence term, which is calculated from differences of velocity components on 
the sublayer mesh. The wall normal velocity V = U 2 is calculated locally at each 
sublayer grid point from differences of the horizontal velocity components using the 
continuity equation, 



The usual boundary layer assumption that dp/dy = 0 is used; hence the pressure 
gradient in (7) for a given horizontal location is set to be constant at all wall-normal 
locations in the sublayer, using the value in the outer flow at the matching point. 
A running time average of the pressure gradient is actually used to smooth the wall 
model. Eq. (7) is discretized and integrated with the same FD2-RK3 scheme used 
in the main code. 

Model JK1. The actual JK model for the inner regions uses velocity scales (w 5 
and Urf) in the eddy viscosity expression (6) that are melds of u T and t/ m , where u m 
is the square root of the maximum Reynolds stress (— tn7) that occurs at a distance 
y max above the wall: 

U 3 = (1 - l)u T + 7 u m , 7 = tanh (y w /£) , £ = u T y mAX /(u r + u m ) , (9a) 

Ud = max(u m ,u r ) . (96) 

Model JK1 calculates Ui from (5), but uses (9) to compute the eddy viscosity in 
(6). In RANS models, u m and y max are determined from the solution of an ODE. 
In LES, the maximum of the stress can in principle be found on the fly at a given 
horizontal position from values of the Reynolds stress in the sublayer and overlying 
outer layer. In practice, this is much more difficult to accomplish with any great 
accuracy, because instantaneous values of the stress along a vertical line fluctuate 
wildly in space and time. Again, a running time average must be used, along with 
some local spatial filtering, in order to smooth the signal to a useful level; then a 
search routine is employed to find the first local maximum of averaged stress moving 
away from the wall at a given horizontal location. Because this is a rather costly 
and cumbersome procedure to employ in LES, its benefits must be shown to be 
substantial to justify its use. 

The computational overheads of the above wall models were about 10, 20, and 
30% of total cost, respectively; however, the number of interior points was halved 
and the time step used was 3 times larger than in a regular, resolved LES, so that 
a savings factor of about 5 was realized. 

The mean streamwise velocities that are obtained using these wall models for 
channel flow are shown in Fig. 2a in comparison with the experimental data (Hussain 
& Reynolds, 1975) and with a LES (Cabot, 1994) for the same parameters with 
the same code without wall models (using 65 wall-normal nodes with about the 
same interior resolution as the LES with wall models and a 64 x 64 horizontal 
mesh). It is seen that there is little difference between the results for different wall 
models in channel flow, suggesting that a simple instantaneous log law provides 


FIGURE 2. Mean streamwise (a) velocity and (b) velocity fluctuation intensi- 
ties in channel flow for LES with different wall models ( JKO, JKOa, 

JK1), compared with a full LES ( , Cabot, 1994) and experimental 

data ( o o o , Hussain & Reynolds, 1975). 
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an adequate, cost-effective wall model in this case. The results for U in the outer 
region are in generally fair agreement with the experimental data and full LES. The 
streamwise velocity fluctuation intensities (u rms ) are shown in Fig. 2b and also show 
fair agreement with experimental and full LES results, with some discrepancies near 
the matching point. Note that there is a large disagreement between the full LES 
results and experimental data in the near-wall region where u rms peaks (y+ < 50). 
The results were insensitive to whether the matching point was at = 32 or 64. 

2. S Boundary layer wall models behind a backward-facing step 

Wall models JKO and JKOa have also been implemented in the LES of flow over 
a backward-facing step using the same FD2-RK3 scheme used for the channel (Ak- 
selvoll & Moin, 1995). The flow has a Reynolds number of 28,000 based on the 
centerline velocity of the inlet flow and the step height h. There is a long inlet 
section 10 h long, 4 h high, and 2 h wide on a 100 x 65 x 96 mesh followed by a 
20 h x bh x 2h outlet section on a 146 x 97 x 96 mesh; both x and y coordinates are 
stretched. The wall model is implemented only along the bottom wall behind the 
step for test purposes, with a 74 x 33 x 48 sublayer mesh embedded below y « 0.073 
or t/ + « 60 at the outlet. No account is taken of the geometry of the corner behind 
the backstep, where there is a weak recirculation zone, but this inaccuracy is not 
expected to affect the bulk of the flow very much. Because only about 10% of the 
grid points are removed from the main calculation and time steps can only be in- 
creased by about 30%, little computational saving is gained from the wall model in 
this case. 

There is a strong adverse pressure gradient between about 3 h and Ih behind 
the step and a concomitant separation bubble in this region. Figure 3 shows the 
near-wall (y/h « 0.10) streamwise pressure gradient from Akselvoll &: Moin’s (1995) 
LES, averaged over time and span; the mean wall-normal gradient of streamwise 
velocity (proportional to the wall stress) is also shown. The assumption that there is 
no wall- normal variation in pressure gradient is found to be good for the most part, 
except in a few regions associated with relatively rapid wall-normal velocities in the 
reattachment region around xjh — 5-8. Preliminary results from the application of 
the JKO wall model (which includes no pressure gradient or advection terms) show 
an underprediction of the level of reversed wall flow (Fig. 3); the recovery region 
around xjh — 10 is also not predicted very well, nor is the recirculation region 
near the step. The level of the post-recovery region near the outlet is predicted 
better; but this region is in fact similar to channel flow or a zero-pressure-gradient 
boundary layer, in which this model was seen to give good results (§2.2). Longer 
runs (currently in progress) are needed to see how the flow adjusts itself further, 
and if the resulting statistically steady flow is predicted adequately. 

The large pressure gradient and advection terms in Eq. (7) are probably required 
to obtain better agreement. For instance, if the streamwise pressure gradient inte- 
grated over the thickness of the sublayer y™, which is about y m dp/dx , is comparable 
to t\ 2 w = vdU jdy\ w , then it can be expected to significantly modify the structure 
of the boundary layer and the wall stress itself. In Fig. 3 the streamwise pressure 
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FIGURE 3. Mean wall-normal gradient of the streamwise velocity at the bottom 

wall dU/dy\ w ( ) and mean streamwise pressure gradient near the bottom 

wall dp /da : ( ) behind a backward-facing step from the LES of Akselvoll & 

Moin (1995). The pressure gradient is scaled by y m /5u, where y m is the thickness of 
the sublayer used in wall model calculations. Preliminary values of mean dU/dy\ w 
predicted with the JKO wall model axe also shown ( • • • ). 

gradient multiplied by y m /5u is seen to be comparable to dU/dy\ w in the separa- 
tion and recirculation regions, and it is likely to have an important effect there. Of 
course, the effect of pressure gradient term will be mitigated to a large degree by 
the advection terms (mostly dU 2 /dx) in the outer part of the sublayer, but these 
terms vanish very near the wall, while the pressure gradient does not. 

Application of the JKOa model, with the addition of large pressure gradient and 
advection terms, shows a much better initial agreement in the reverse flow region, 
although the recovery region around x/h — 10 is still not well predicted. The 
region around x/h = 5 near the head of the separation bubble in the reattachment 
zone, characterized by downflows that axe strong in comparison with horizontal 
flow, has led to numerical instability in the sublayer calculation. The cause of 
this is still not known, but it appears to be associated with very large advection 
terms d(UiV)/dy at locations of rapid downflow. These are also regions where the 
assumption of constant horizontal pressure gradients breaks down and the boundary 
layer equations are known to be invalid. 

3. Future plans 

Some fundamental tests need to be performed on backward-facing step flow fields 
near the bottom wall, such as the momentum balance at different scales that was 
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performed for channel flow (§2.1). DNS and LES fields will be studied to attempt 
to determine, for example, how the changes in pressure gradient affect in detail the 
wall stress and what terms in the momentum equation are most important in the 
regions of strong downflow at the head of the separation bubble. 

LES with the simple JKO wall model (essentially the smooth meld of the log 
law and viscous law) will be run over long times to statistical equilibrium to get 
a fair assessment of that model’s performance. The same model with advection 
and running time-averaged pressure terms (JKOa) will also be run to longer times 
if the present numerical instability can be cured. An attempt will also be made to 
implement the JK1 wall model in the backward-facing step flow, which requires a 
determination of the maximal shear stress (averaged in some sense) above the wall in 
order to determine a model velocity scale. Search routines like that used in channel 
flow, and perhaps a curve fitting scheme applied to the shear stress profiles, will be 
tried; however, there is always some arbitrariness in these approaches. Alternative, 
more easily determined, and better quantified velocity scales will also be considered. 
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